Soliton creation during a Bose-Einstein condensation 
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We use stochastic Gross-Pitaevskii equation to study dynamics of Bose-Einstein condensation. 
We show that coohng into a Bose-Einstein condensate (BEG) can create sohtons with density given 
by the cooUng rate and by the critical exponents of the transition. Thus, counting soHtons left in 
its wake should allow one to determine the critical exponents z and v for a BEG phase transition. 
The same information can be extracted from two-point correlation functions. 



When temperature of a bosonic cloud is lowered, atoms 
undergo a thermal gas - Bose-Einstein condensate phase 
transition. Sufficiently slow cooling process creates a con- 
densate in the ground state. Its macroscopic quantum 
properties have been studied over the last decade. Faster 
cooling should result in creation of an excited condensate 
that contains vortices and solitons. 

Non-equilibrium condensation is an example of dynam- 
ics of phase transitions. Their study started with in- 
vestigations of topological defect creation in cosmologi- 
cal phase transitions fl]. It was then pointed out that 
condensed matter systems (superconductors, liquid "^He 
and ''He, etc.) can be used to test this cosmological sce- 
nario, and that the nature of the transition will determine 
density of defects [2]. Kibble-Zurek mechanism (KZM) 
describes how topological defects (vortices, monopoles, 
kinks, etc.) are created during non-equilibrium second 
order phase transitions (see for experimental evi- 

dence) . 

Data on spontaneous creation of solitons and vortices 
during Bose-Einstein condensation are still scarce. Two 
groups have seen them: the Anderson group observed 
vortices in a three dimensional trap [3], while the En- 
gels' group reported observation of solitons in a quasi 
one dimensional (ID) setup that does not support topo- 
logical defects Q. Therefore, it is important to under- 
stand non-equilibrium dynamics of Bose-Einstein con- 
densation. Moreover, it is also interesting to find out 
if the KZM, studied theoretically and experimentally 
in various physical systems 0-11, applies to non- 

topological excitations (e.g., solitons). We propose to 
kill two birds with one stone by studying dynamics of 
Bose-Einstein condensation in a quasi-one ID system. 

As the atom cloud is cooled from above the critical 
point, the system adjusts adiabatically to driving because 
its relaxation time is initially short. Near the critical 
point, however, the relaxation time diverges due to criti- 
cal slowing down. The system goes out of equilibrium be- 
fore reaching the critical point and approximately freezes 
out entering the impulse stage of its dynamics. As the 
condensate forms, its different parts choose to break the 
symmetry in an uncorrelated way, which results in cre- 
ation of defects. Their spacing depends on the quench 
rate given by the rate of temperature change [2|: the 
slower we go, the more adiabatic the evolution is, and 



so the larger the separation between the defects will be. 
The transition from adiabatic to impulse regime happens 
when the relaxation time r becomes comparable to the 
quench rate e/i {e is the distance from the critical point): 



(1) 



Assuming a linear quench, i.e., \i\ = Tq^ , where tq is the 
quench timescale, the system goes out of equilibrium at 



(2) 



This imprints a characteristic time scale t = Tgi ~ 
Tq^'^^^^^^ onto the system, so that defect density de- 
pends on time through tjt or equivalently e{t)/e. The 
typical distance between defects scales as 



Q 



(3) 

in the KZM freeze out picture Above z and v are 
the critical exponents defining the equilibrium coherence 
length ^ and relaxation timescale r: 



(4) 



Thus, KZM (presented above) proposes a way to extract 
essential features of the non-equilibrium phase transition 
dynamics from equilibrium critical behavior of a system. 

The exact simulation of Bose-Einstein condensation is 
extremely involved (if at all possible) [l^]- We restrict 
ourselves to a tractable model: the stochastic Gross- 
Pitaevskii equation (SGPE) that has been successfully 
applied to studies of Bose-Einstein condensates lately 
[lO|. Our calculations provide the first numerical results 
on dynamics of soliton production in the course of second 
order phase transitions, and extend the KZM to non- 
topological excitations. Former studies of ID models, 
driven across a critical point, were focused on systems 
supporting topological defects (e.g., kinks [1]). 



Stochastic Gross-Pitaevskii equation reads : 



{i - j)dt(f> ■ 



-dl^ + ecf, + g\<l>\^^ + ^{x,t), 



where noise, coming from a thermal cloud, satisfies 
{d{x, t)r {x', t')) = 2-fTS{x - x')6{t - t'). 



(5) 



(6) 
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Above 7 represents damping coming from thermal cloud 
- condensate interactions, T is the thermal cloud tem- 
perature, and £ = — /i, where [i is the chemical potential. 

Forgetting for a while about damping and noise we see 
that the system is described by the energy functional 



E = 



dx - \dx 



-Vm), VM)^e\cj^f + ^\cb\\ (7) 



For £ > (/I < 0) the minimum of ^^df/il) corresponds 
to (j) — 0: order parameter vanishes above condensation 
temperature, and the system is in the symmetric phase. 
When £ < (/i > 0) there is a minimum of ^(|(/>|) at 
|0p — — |, and the system is in the broken-symmetry 

phase where <p = v^Tff exp(i0): a condensate forms. 

We quench the system from the symmetric to the 
broken-symmetry phase. The phase 9 will not be cho- 
sen uniformly in the broken-symmetry phase because the 
quench freezes a characteristic coherence length, ^, when 
the system goes out of equilibrium on the symmetric side. 
This leads to phase gradients that seed soliton-like den- 
sity notches around which the phase of the order param- 
eter changes abruptly [l2. 13|. 

For simplicity, we model cooling of a bosonic gas by 
assuming that chemical potential grows across the criti- 
cal point, while damping 7 and noise correlations remain 
constant. This approximation is motivated by the fact 
that defect production takes place near the critical point, 
where chemical potential changes drive the transition in- 
fluencing dynamics most, while the variations in 7 and 
noise correlators are of secondary importance. 

In the following we need freeze out values of i and ^. 
As z = 2 and = 1/2 for SGPE, Eqs. © and @ imply 



//4 



(8) 



We drive the system from the symmetric phase to the 
broken-symmetry phase at a constant rate: 



eit) = 



t 



(9) 



where t = —Sqtq — )■ SqTq. The evolution starts away 
from the critical point at £ = £0 3> 1. The system passes 
through the critical point at e,t — 0, and ends its evolu- 
tion far away from the critical point at £ = — £0. 

The discussion of KZM can be analytically illustrated 
in our model on the symmetric side. After dropping the 
nonlinear term we are left with Langevin-type equation: 



{i - -f)dt(j) ■ 



1 



dl(t>'Ve(t> + 'd{x,t). 



(10) 



This equation can be solved exactly providing the follow- 
ing relaxation time and coherence length in equilibrium: 



1 + 7^1 , 

T ^ , C 

7 £ 



1 1 



(11) 



2 




L (a) 








direction of evolution 








> 






- , 1 , 1 , 1 , - 




100 



FIG. 1: Half-width of the averaged two-point correlation func- 
tion (|15p during the quench: ^1/2. (a) and (b): ^1/2 before 
and after rescalings. On both plots black, red, and green 
lines correspond to tq — 5, 25 and 85, respectively. The in- 
set magnifies the region where the KZM works best, (c) Il- 
lustration that quench imprints a characteristic length scale 



'T'o^ We study ^1/2 at the critical point. 



The scale is logarithmic on both axes: circles show numer- 
ics, while the straight line is a fit to numerics: ln^i/2(e = 
0) = -1.53 ± 0.03 + (0.24 ± 0.01) Inrg. The fit confirms 
the predicted 1/4 exponent, (d) Illustration that quench im- 
prints a characteristic time scale i ~ ^^"'/{i+z") _ rg'^. We 
study time At passing from entering the broken-symmetry 
phase to the moment when ^1/2 /£, reaches a threshold, as- 
sumed here to be 0.5, i.e., to be near the center of the re- 
gion where rescalings work on the broken-symmetry side [see 
inset of Fig. (b)]. The scale is logarithmic on both axes: 
circles show numerics, while straight line is a fit to numerics: 
In At = 2.60 ±0.04 +(0.49 ±0.01) Inrg. This is in good agree- 
ment with theoretical prediction of At ~ '''q^ ■ We use for this 

plot 5 taken from the fit done in plot (c): ^ — tq^^. Fits from 
plots (c) and (d) allow for measurement of z and 1/ critical 
exponents. All results presented above are averages over 100 
runs with different noise. See [Ij for other parameters. 



respectively. It implies z = 2 and ly — 1/2. We obtained 
pT|) by studying the two point correlation function 

Cix,t\x',t') = {cj,{x,t)rix',t')) - {^{x,t)){rix',t')), 

(12) 

where (...) denotes averaging over different noise real- 
izations. We found that C{x, t\x' , t) decays on the length 
scale ^ while C(a:, t|a;, t') decays on the timescale t. 

Solving ID) with (0) and PT|) we determine that the 
adiabatic - impulse border is at 



£ = 



'l±72 1 
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so the system "freezes" when its coherence length equals 

Assuming for simplicity that Sq ^ oo, and solving an- 
alytically ^ with ^ we find that far away from the 
critical point the system follows adiabatically the instan- 
taneous equilibrium solution 



C{x,t\x\t) = 



T 



exp 



[-\x^x'\VW))- (13) 



Near the critical point, i.e., when e < e we have to 
refer to the exact solution that reads 

C{x,t\x',t) = {\m\')ej{\x~x'\/ieit)/e), (14) 

where /(a, 6) — dfc cos(fca) cxp((6+fc^)^)erfc(6-|-fc^), 
(■■)eq denotes equilibrium averaging, and (|0(e)p)eg = 
T/\/2i. Naturally, ^ reduces to (HH) for e » e. 

The result shows that correlations are induced 
on the length scale which is given by the correla- 
tion length at the border between adiabatic and impulse 
regimes. Another interesting feature of ([T4l) is that it 
scales as Tq"^: the equilibrium correlations calcu- 
lated at e(t) = e scale in the same way. Similarly, there is 
a freeze out time scale t imprinted into system's dynam- 
ics through the relation e{t)/e — t/t. All these results 
are in perfect agreement with KZM. 

On the broken-symmetry side we rely on numerics. As 
above, we test our theory by looking at the characteris- 
tic length and time scales imprinted onto the system by 
the quench. This is done by studying half-width of the 
averaged two-point correlation function, i.e.. 



1 



dx \C{x, t\x + r,t)\ , 



(15) 



where the averaging over the system size I is applied. Fig. 
[1] shows that in whole symmetric phase, as well as in the 

1/2 

broken-symmetry phase for —stq < 25, the following 
scaling works very well 



1/2 



^1/2 -^f{e/i) 

(16) 

Scaling (|16p is in perfect agreement with KZM and 
allows for measurement of critical exponents z and h' 
(see Fig. [T]). We use here SGPE mean- field exponents 
{z — 2 and v — 1/2) to compare our numerics to the- 
ory. In an actual experiment, however, other values of 
critical exponents may be relevant (see 15| and refer- 
ences cited therein). Indeed, critical exponents can be 
experimentally studied, for example, with the help of 
cavity- assisted cold atom counting recently explored in 
the Esslinger's group at ETH Zurich |15|]. There, the 




FIG. 2: Snapshot of density |(;/>(a::)|^, phase arg(<jf>), and gra- 
dient of phase ^arg(<^) (the velocity field). The snapshot 
is taken for one evolution (no averaging) at e{t) = 10 for 
TQ — 10. See for other parameters. 



system was near equilibrium and so the study of two- 
point correlation functions was aimed at determination 
of the critical exponent v. Nonequilibrium version of this 
experiment can provide a new way for measurement of 
the exponent v and can reveal the dynamical exponent 
z. All this is possible due to presence of the character- 
istic length scale ^1/2 in the non-equilibrium state of the 
condensate. The most striking consequence of existence 
of this length scale is the creation of solitons. 

Qualitatively, we observe soliton-likc solutions in the 
broken-symmetry phase: see Fig. [2] There are several 
deep density notches there, and the phase of the order pa- 
rameter changes steeply around them. These are typical 
signatures of solitons [1^ . They are in qualitative agree- 
ment with Engels' experiments studying density profiles 
after the condensation process [5] . 

Quantitatively, we would like to find out if the criti- 
cal scalings ^ can be retrieved from soliton counting. 
This is of both fundamental and practical interest: soli- 
ton counting can be significantly easier than measure- 
ment of the correlation functions. The typical number 
of solitons shall be inversely proportional to the size of 
correlated domains that have chosen to break the symme- 
try in the same way: phase jumps between the domains 
provide seeds for solitons. Thus, we predict 



# of solitons ^ ^-^ 



/2 



/2- 'r"f[err")--Q'''fH 

(17) 

We count solitons by fitting the solitonic solution of the 
Gross-Pitaevskii equation (no noise/damping) around 
every minimum of the density |0p: 

+ {no ) tanh^[(a; - XQ)^/g(j^^~^n^^^^]. 
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of the interacting Bose-Einstein condensate. 

Summarizing, we have shown that two-point correla- 
tion functions after non-equilibrium quench encode crit- 
ical exponents of the system. We have also predicted 
that non-equilibrium condensation results in creation 
of soliton-like excitations whose number depends on a 
quench rate via the critical exponents. Our results can 
be "inverted" and used for the experimental determina- 
tion of the critical exponents for the normal gas - Bose- 
Einstein condensate phase transition. 

We are grateful to Peter Engels for showing us his un- 
published experimental data and for stimulating discus- 
sions. We thank Ashton Bradley for his very useful com- 
ments. We acknowledge the support of the U.S. Depart- 
ment of Energy through the LANL/LDRD Program. 



FIG. 3: Number of solitons in the broken-symmetry phase 
before and after rescalings. On both plots black, red and 
green lines correspond to tq — 10, 40 and 70, respectively. 
The inset magnifies the region where KZM works best. The 
results present average over 100 runs with different noise. See 
[3| for other parameters. 



We use only modulus of the order parameter as it is typ- 
ically the only quantity measurable in a standard BEG 
setup. The fit gives the three parameters of the soliton so- 
lution: position (xq), minimum density (rimm), and back- 
ground density (no)- The latter two are related to soliton 
and sound velocities, respectively. We then compare the 
fitted depth of the soliton, uq — rimm, to the numerical 
data. When the two do not differ by more than 50%, 
we count the density minimum as a soliton. We have 
checked that the same conclusions are obtained for other 
reasonable thresholds between 30% and 60%. The out- 
come of this procedure is presented in Fig. [3l where the 
number of solitons for —£Tq < 25 follows closely (fT7|) 
in agreement with KZM and our studies of the width of 
the two-point correlation functions. The drawback with 
respect to the latter is that solitons need some time to 
develop before they can be counted, thus their macro- 
scopic number can be observed in a narrow window of 

1/2 

20 < ~sTq < 25. Another complication is that solitons 
are more vulnerable to damping than two-point correla- 
tion function: damping removes solitons from the system 
by decreasing their depth. This shall not be a problem 
with an actual experiment as there is negligible damp- 
ing/noise by the end of condensation 

Our results were obtained for atoms in a quasi-lD ho- 
mogeneous configuration that can be experimentally re- 
alized 16|. In harmonic traps with cigar-shaped BEC's 
homogeneous scalings would be modified by causality- 
related considerations 0. In either case, the experi- 
mental data will shed light on the real critical exponents 
of the cold atom system and facilitate the first experi- 
mental determination of the dynamical, i.e., z exponent 
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